Reuse solver residual in calibrator - #231
Conversation
The block update_gradient methods previously hand-coded the residual r = E*ydot + F*y + c, duplicating the residual the solver already assembles in SparseSystem::update_residual (r = -C - E*ydot - F*y). Reuse the solver residual during calibration and have update_gradient assemble only the Jacobian of the residual with respect to the parameters. The Jacobian signs are flipped to match the solver's sign convention; the overall sign cancels in the LM normal equations, so the optimization is unchanged. - LevenbergMarquardtOptimizer now sets the model parameters to the current alpha, assembles E/F once, and computes the residual per observation via update_solution + update_residual. - calibrate registers the parameters with the model so the solver's assembly can read their values. - Blocks without calibration parameters (e.g. a normal junction) are skipped in the Jacobian loop, so Junction::update_gradient is removed. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ndeleev-a7e408 # Conflicts: # src/model/Junction.h # src/optimize/LevenbergMarquardtOptimizer.cpp # src/optimize/calibrate.cpp
mrp089
changed the title
|
@federicaninno, can you review this? |
|
@mrp089 sure, will do! |
|
|
||
| // Set up the solver system so that its residual assembly can be reused for a | ||
| // single observation at a time. | ||
| system = SparseSystem(num_vars); |
There was a problem hiding this comment.
| system = SparseSystem(num_vars); | |
| system = SparseSystem(std::max(num_vars, num_eqns)); |
The calibrator's model omits boundary condition blocks (their values come from the observations instead), so num_eqns and num_vars need not match here as they do for a full simulation. Size the system by the larger of the two so that both equation-row and variable-column indices fit.
| // here. | ||
| model->update_solution(system, y_dpoint, dy_dpoint); | ||
| system.update_residual(y_dpoint, dy_dpoint); | ||
| residual.segment(num_eqns * i, num_eqns) = system.residual; |
There was a problem hiding this comment.
| residual.segment(num_eqns * i, num_eqns) = system.residual; | |
| residual.segment(num_eqns * i, num_eqns) = system.residual.head(num_eqns); |
system.residual has max(num_vars, num_eqns) slots (see comment at line 32), which is num_vars in the calibrator, since there are no boundary condition blocks. However, only the first num_eqns slots ever get filled in, the rest stay zero. .head(num_eqns) keeps just those first num_eqns slots so the result fits into the num_eqns-sized segment.
federicaninno
left a comment
federicaninno
left a comment
There was a problem hiding this comment.
Hi @mrp089, I tried running the calibrator on this branch against the existing test case (steadyFlow_calibration.json), and I got this error:
"Assertion failed: (rows == this->rows() && cols == this->cols() && "DenseBase::resize() does not actually allow to resize."), function resize, file DenseBase.h, line 262.
zsh: abort ./svzerodcalibrator ../tests/cases/steadyFlow_calibration.json try.json"
It seems like "calibrate.cpp" only builds vessel/junction blocks and skips the boundary-condition blocks (FlowReferenceBC, ResistanceBC, etc.), since those values come straight from the observation data instead of being solved for. But those BC blocks are normally what add the extra equations that make the system square (num_vars == num_eqns). Without them, we end up with more variables than equations, which then creates problems in "LevenbergMarquardtOptimizer.cpp". At line 32: system = SparseSystem(num_vars); sizes system.residual to num_vars. Later, at line 102: residual.segment(num_eqns * i, num_eqns) = system.residual tries to stuff that num_vars-sized vector into a num_eqns-sized slot, and cannot resize to fit, hence the crash. I am leaving a possible fix directly as comments on the two lines in question, which I tested and worked.
Let me know what you think!
2 participants
Current situation
Closes #153.
Release Notes
The calibrator reuses the code from the solver to calculate the residual, which previously needed to be reimplemented for each block. This makes the code more straightforward and makes calibrating new blocks easier.
Documentation
No change in functionality
Testing
Should work as before
Code of Conduct & Contributing Guidelines